Survival MANOVA
Kevin Middleton
3/7/2017
Setup
set.seed(234876)
library(MCMCglmm)## Loading required package: Matrix## Loading required package: coda## Loading required package: apelibrary(tidyverse)## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr## Conflicts with tidy packages ----------------------------------------------## expand(): tidyr, Matrix
## filter(): dplyr, stats
## lag(): dplyr, statslibrary(parallel)
library(knitr)
source("h2life_load_data.R")## Parsed with column specification:
## cols(
## setDate = col_date(format = ""),
## flipDate = col_date(format = ""),
## days = col_integer(),
## age = col_integer(),
## NewAge = col_integer(),
## fID = col_character(),
## id = col_character(),
## sireid = col_character(),
## damid = col_character(),
## repl = col_character(),
## treat = col_character(),
## NstartF = col_integer(),
## status = col_integer()
## )genet_corr <- tibble(model = character(),
HS_STD = numeric(),
LY_STD = numeric(),
HS_LY = numeric(),
n_eff = numeric())
iter <- 2.5e5
burnin <- 2e4
thinning <- 500
chains <- 12
rerun <- TRUEMANOVA analysis
Poster for fly meeting
if (rerun) {
HS <- h2life %>%
filter(treat == "HS") %>% rename(Age_HS = NewAge) %>%
as.data.frame()
LY <- h2life %>%
filter(treat == "LY") %>% rename(Age_LY = NewAge) %>%
as.data.frame()
STD <- h2life %>%
filter(treat == "STD") %>% rename(Age_STD = NewAge) %>%
as.data.frame()
h2life_mrg <- full_join(full_join(HS, LY), STD)
prior1 <- list(R = list(V = diag(3) * 1.002, nu = 1.002),
G = list(G1 = list(V = diag(3) * 1.002, nu = 0.002)))
prior2 <- list(R = list(V = diag(3) / 3, nu = 1.002),
G = list(G1 = list(V = diag(3) / 3, nu = 0.002)))
prior3 <- list(R = list(V = diag(3) * 1.002, nu = 1),
G = list(G1 = list(V = diag(3) * 0.5, nu = 0.002)))
priors <- list(prior1, prior2, prior3)
prior_names <- c("Tri: V = diag(3) * 1.002, nu = 0.02",
"Tri: V = diag(3) / 3, nu = 0.02",
"Tri: V = diag(3) * 0.5, nu = 0.02")
for (ii in 1:length(priors)) {
prior <- priors[[ii]]
fm <- mclapply(1:chains, function(i) {
MCMCglmm(cbind(Age_STD, Age_HS, Age_LY) ~ trait - 1,
random = ~ us(trait):animal,
rcov = ~ idh(trait):units,
data = h2life_mrg,
prior = prior,
pedigree = pedigree,
family = rep("gaussian", 3),
nitt = iter,
burnin = burnin,
thin = thinning)
}, mc.cores = chains)
outfile <- paste0("tri_model_prior", ii, ".Rda")
save(fm, file = outfile)
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
re <- as.mcmc(as.matrix(re))
n_eff <- effectiveSize(re)
# STD vs. HS
HS_STD <- re[ , "traitAge_HS:traitAge_STD.animal"] /
sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
re[ , "traitAge_HS:traitAge_HS.animal"])
# STD vs. LY
LY_STD <- re[ , "traitAge_LY:traitAge_STD.animal"] /
sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
re[ , "traitAge_LY:traitAge_LY.animal"])
# LY vs. HS
HS_LY <- re[ , "traitAge_HS:traitAge_LY.animal"] /
sqrt(re[ , "traitAge_LY:traitAge_LY.animal"] *
re[ , "traitAge_HS:traitAge_HS.animal"])
genet_corr <- bind_rows(genet_corr,
tibble(model = prior_names[[ii]],
HS_STD = median(HS_STD),
LY_STD = median(LY_STD),
HS_LY = median(HS_LY),
n_eff = mean(n_eff)))
}
write_csv(genet_corr, path = "Genetic_Correlations.csv")
}## Joining, by = c("sire", "dam", "animal", "treat")
## Joining, by = c("sire", "dam", "animal", "treat")Analyze model
load("tri_model_prior1.Rda")
fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe, ask = FALSE)
plot(fe[, 1, drop = FALSE], ask = FALSE)
plot(fe[, 2, drop = FALSE], ask = FALSE)
plot(fe[, 3, drop = FALSE], ask = FALSE)
# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
plot(re[, 1, drop = FALSE], ask = FALSE)
plot(re[, 2, drop = FALSE], ask = FALSE)
plot(re[, 3, drop = FALSE], ask = FALSE)
autocorr.diag(re)## traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## Lag 0 1.0000000000 1.000000000
## Lag 500 0.0842692033 0.076989902
## Lag 2500 -0.0178249648 -0.018361472
## Lag 5000 0.0151468155 -0.001760693
## Lag 25000 0.0008825723 -0.007553334
## traitAge_LY:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## Lag 0 1.000000000 1.000000000
## Lag 500 0.069088642 0.076989902
## Lag 2500 -0.024071358 -0.018361472
## Lag 5000 0.009317922 -0.001760693
## Lag 25000 0.003487370 -0.007553334
## traitAge_HS:traitAge_HS.animal traitAge_LY:traitAge_HS.animal
## Lag 0 1.000000000 1.000000000
## Lag 500 0.083773312 0.083790353
## Lag 2500 -0.004399008 0.004611076
## Lag 5000 -0.022509111 -0.008720982
## Lag 25000 -0.016121417 -0.007318574
## traitAge_STD:traitAge_LY.animal traitAge_HS:traitAge_LY.animal
## Lag 0 1.000000000 1.000000000
## Lag 500 0.069088642 0.083790353
## Lag 2500 -0.024071358 0.004611076
## Lag 5000 0.009317922 -0.008720982
## Lag 25000 0.003487370 -0.007318574
## traitAge_LY:traitAge_LY.animal traitAge_STD.units
## Lag 0 1.000000000 1.000000000
## Lag 500 0.083268922 0.075959690
## Lag 2500 0.005213926 -0.012102331
## Lag 5000 0.007473910 0.015920811
## Lag 25000 -0.014482780 0.004354346
## traitAge_HS.units traitAge_LY.units
## Lag 0 1.000000000 1.0000000000
## Lag 500 0.075349367 0.0642032980
## Lag 2500 0.007898546 -0.0153543256
## Lag 5000 -0.024638971 -0.0001941856
## Lag 25000 -0.019878975 -0.0175330898effectiveSize(re)## traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## 5004.529 5119.084
## traitAge_LY:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## 5161.457 5119.084
## traitAge_HS:traitAge_HS.animal traitAge_LY:traitAge_HS.animal
## 4858.866 4784.512
## traitAge_STD:traitAge_LY.animal traitAge_HS:traitAge_LY.animal
## 5161.457 4784.512
## traitAge_LY:traitAge_LY.animal traitAge_STD.units
## 4899.712 4807.313
## traitAge_HS.units traitAge_LY.units
## 4849.351 5019.077# Concatenate samples
re <- as.mcmc(as.matrix(re))
head(re)## Markov Chain Monte Carlo (MCMC) output:
## Start = 1
## End = 7
## Thinning interval = 1
## traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## [1,] 0.4146686 0.12177726
## [2,] 0.4518303 0.11549431
## [3,] 0.3258092 0.09183121
## [4,] 0.3272157 0.05843721
## [5,] 0.6759533 0.26814659
## [6,] 0.5333801 0.23382353
## [7,] 0.3621810 0.15351107
## traitAge_LY:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## [1,] 0.10732055 0.12177726
## [2,] 0.05899287 0.11549431
## [3,] 0.18131683 0.09183121
## [4,] 0.15226261 0.05843721
## [5,] 0.25923669 0.26814659
## [6,] 0.24500887 0.23382353
## [7,] 0.15202631 0.15351107
## traitAge_HS:traitAge_HS.animal traitAge_LY:traitAge_HS.animal
## [1,] 0.3561362 0.10024041
## [2,] 0.2730102 0.05165598
## [3,] 0.3557770 0.13505231
## [4,] 0.1921802 0.09424014
## [5,] 0.3621698 0.10099553
## [6,] 0.4371556 0.30262471
## [7,] 0.4382596 0.24309716
## traitAge_STD:traitAge_LY.animal traitAge_HS:traitAge_LY.animal
## [1,] 0.10732055 0.10024041
## [2,] 0.05899287 0.05165598
## [3,] 0.18131683 0.13505231
## [4,] 0.15226261 0.09424014
## [5,] 0.25923669 0.10099553
## [6,] 0.24500887 0.30262471
## [7,] 0.15202631 0.24309716
## traitAge_LY:traitAge_LY.animal traitAge_STD.units traitAge_HS.units
## [1,] 0.1728861 0.5121829 0.5017728
## [2,] 0.2074450 0.5128448 0.5476254
## [3,] 0.3866927 0.5540212 0.5016486
## [4,] 0.2524715 0.5582436 0.5875449
## [5,] 0.2146196 0.3887915 0.4857655
## [6,] 0.3120685 0.4636957 0.4390513
## [7,] 0.2938009 0.5428062 0.4631711
## traitAge_LY.units
## [1,] 0.6115914
## [2,] 0.5809705
## [3,] 0.4701095
## [4,] 0.5232854
## [5,] 0.5458642
## [6,] 0.5432736
## [7,] 0.5280122# STD vs. HS
plot(re[ , "traitAge_HS:traitAge_STD.animal"])
plot(re[ , "traitAge_STD:traitAge_STD.animal"])
HS_STD <- re[ , "traitAge_HS:traitAge_STD.animal"] /
sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
re[ , "traitAge_HS:traitAge_HS.animal"])
plot(HS_STD)
median(HS_STD)## [1] 0.2917385HPDinterval(HS_STD)## lower upper
## var1 0.006316709 0.5387122
## attr(,"Probability")
## [1] 0.95# STD vs. LY
LY_STD <- re[ , "traitAge_LY:traitAge_STD.animal"] /
sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
re[ , "traitAge_LY:traitAge_LY.animal"])
plot(LY_STD)
median(LY_STD)## [1] 0.4928118HPDinterval(LY_STD)## lower upper
## var1 0.2525297 0.7072765
## attr(,"Probability")
## [1] 0.95# LY vs. HS
HS_LY <- re[ , "traitAge_HS:traitAge_LY.animal"] /
sqrt(re[ , "traitAge_LY:traitAge_LY.animal"] *
re[ , "traitAge_HS:traitAge_HS.animal"])
plot(HS_LY)
median(HS_LY)## [1] 0.348081HPDinterval(HS_LY)## lower upper
## var1 0.06371274 0.5937976
## attr(,"Probability")
## [1] 0.95Plot for poster
library(tidyverse)
library(cowplot)
M <- data_frame(`HS vs. STD` = as.numeric(HS_STD),
`LY vs. STD` = as.numeric(LY_STD),
`HS vs. LY` = as.numeric(HS_LY))
colMeans(M)## HS vs. STD LY vs. STD HS vs. LY
## 0.2864158 0.4844311 0.3401840M %>% gather(Comparison, value) %>%
ggplot(aes(value, color = Comparison)) +
geom_line(stat = "density", size = 2) +
labs(x = "Genetic Correlation", y = "Density") +
theme(legend.position = c(0.15, 0.85),
text = element_text(size = 24),
legend.text = element_text(size = 18))
ggsave(last_plot(), file = "Genetic_Correlation_Plot.pdf",
width = 8, height = 6)Model following Ingleby et al. 2013
Setup data as above
V_nu_to_a_b <- function(V, nu) {
require(tidyverse)
require(invgamma)
a <- nu / 2
b <- (nu * V) / 2
p <- tibble(
x = seq(0, 10, length = 200),
y = dinvgamma(x, a, b)) %>%
ggplot(aes(x, y)) +
geom_line()
print(p)
return(list(alpha = a, beta = b))
}V_nu_to_a_b(1, 0.002)## Loading required package: invgamma## Warning: Removed 1 rows containing missing values (geom_path).
## $alpha
## [1] 0.001
##
## $beta
## [1] 0.001genet_corr <- read_csv("Genetic_Correlations.csv")## Parsed with column specification:
## cols(
## model = col_character(),
## HS_STD = col_double(),
## LY_STD = col_double(),
## HS_LY = col_double(),
## n_eff = col_double()
## )if (rerun) {
prior1 <- list(R = list(V = 1, nu = 0.002),
G = list(G1 = list(V = diag(3) / 3, nu = 0.002)))
prior2 <- list(R = list(V = 1, nu = 0.002),
G = list(G1 = list(V = diag(3) * 0.5 * 0.7, nu = 0.002)))
M <- matrix(0.5 * 0.7, 3, 3)
diag(M) <- 0.5
prior3 <- list(R = list(V = 1, nu = 0.002),
G = list(G1 = list(V = M, nu = 0.002)))
priors <- list(prior1, prior2, prior3)
prior_names <- c("Sire: V = diag(3) / 3, nu = 0.002",
"Sire: V = diag(3) * 0.5 * 0.7, nu = 0.002",
"Sire: V = 0.5 diag, 0.35 off-diag, nu = 0.002")
for (ii in 1:length(priors)) {
prior <- priors[[ii]]
fm <- mclapply(1:chains, function(i) {
MCMCglmm(NewAge ~ treat,
random = ~ us(treat):sire,
data = h2life,
prior = prior,
family = "gaussian",
nitt = iter,
burnin = burnin,
thin = thinning,
verbose = TRUE)
}, mc.cores = chains)
outfile <- paste0("sire_model_prior", ii, ".Rda")
save(fm, file = outfile)
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
re <- as.mcmc(as.matrix(re))
n_eff <- effectiveSize(re)
# STD vs. HS
HS_STD <- re[ , "treatHS:treatSTD.sire"] /
sqrt(re[ , "treatHS:treatHS.sire"] *
re[ , "treatSTD:treatSTD.sire"])
median(HS_STD)
# STD vs. LY
LY_STD <- re[ , "treatSTD:treatLY.sire"] /
sqrt(re[ , "treatSTD:treatSTD.sire"] *
re[ , "treatLY:treatLY.sire"])
median(LY_STD)
# LY vs. HS
HS_LY <- re[ , "treatHS:treatLY.sire"] /
sqrt(re[ , "treatLY:treatLY.sire"] *
re[ , "treatHS:treatHS.sire"])
median(HS_LY)
genet_corr <- bind_rows(genet_corr,
tibble(model = prior_names[[ii]],
HS_STD = median(HS_STD),
LY_STD = median(LY_STD),
HS_LY = median(HS_LY),
n_eff = mean(n_eff)))
}
write_csv(genet_corr, path = "Genetic_Correlations.csv")
kable(genet_corr, digits = 3)
}| model | HS_STD | LY_STD | HS_LY | n_eff |
|---|---|---|---|---|
| Tri: V = diag(3) * 1.002, nu = 0.02 | 0.298 | 0.494 | 0.347 | 4776.086 |
| Tri: V = diag(3) / 3, nu = 0.02 | 0.295 | 0.499 | 0.352 | 4820.542 |
| Tri: V = diag(3) * 0.5, nu = 0.02 | 0.294 | 0.491 | 0.351 | 4812.753 |
| Sire: V = diag(3) / 3, nu = 0.002 | 0.515 | 0.486 | 0.562 | 5436.460 |
| Sire: V = diag(3) * 0.5 * 0.7, nu = 0.002 | 0.514 | 0.493 | 0.566 | 5499.585 |
| Sire: V = 0.5 diag, 0.35 off-diag, nu = 0.002 | 0.520 | 0.489 | 0.566 | 5520.000 |
Analyze model
load("sire_model_prior1.Rda")
fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe, ask = FALSE)
# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)
plot(re[, 1:3], ask = FALSE)
plot(re[, 4:6], ask = FALSE)
plot(re[, 7:9], ask = FALSE)
plot(re[, 10], ask = FALSE)
effectiveSize(re)## treatHS:treatHS.sire treatLY:treatHS.sire treatSTD:treatHS.sire
## 5308.632 5520.000 5565.375
## treatHS:treatLY.sire treatLY:treatLY.sire treatSTD:treatLY.sire
## 5520.000 5274.867 5358.924
## treatHS:treatSTD.sire treatLY:treatSTD.sire treatSTD:treatSTD.sire
## 5565.375 5358.924 5534.616
## units
## 5358.443# Concatenate samples
re <- as.mcmc(as.matrix(re))
colnames(re)## [1] "treatHS:treatHS.sire" "treatLY:treatHS.sire"
## [3] "treatSTD:treatHS.sire" "treatHS:treatLY.sire"
## [5] "treatLY:treatLY.sire" "treatSTD:treatLY.sire"
## [7] "treatHS:treatSTD.sire" "treatLY:treatSTD.sire"
## [9] "treatSTD:treatSTD.sire" "units"# ??? Heritability ????
HS <- re[, "treatHS:treatHS.sire"] / (re[, "treatHS:treatHS.sire"] + re[, "units"])
median(HS)## [1] 0.1500825LY <- re[, "treatLY:treatLY.sire"] / (re[, "treatLY:treatLY.sire"] + re[, "units"])
median(LY)## [1] 0.09315043STD <- re[, "treatSTD:treatSTD.sire"] / (re[, "treatSTD:treatSTD.sire"] + re[, "units"])
median(STD)## [1] 0.1680397# STD vs. HS
HS_STD <- re[ , "treatHS:treatSTD.sire"] /
sqrt(re[ , "treatHS:treatHS.sire"] *
re[ , "treatSTD:treatSTD.sire"])
plot(HS_STD)
median(HS_STD)## [1] 0.5150163HPDinterval(HS_STD)## lower upper
## var1 0.1742536 0.8009837
## attr(,"Probability")
## [1] 0.95# STD vs. LY
LY_STD <- re[ , "treatSTD:treatLY.sire"] /
sqrt(re[ , "treatSTD:treatSTD.sire"] *
re[ , "treatLY:treatLY.sire"])
plot(LY_STD)
median(LY_STD)## [1] 0.486466HPDinterval(LY_STD)## lower upper
## var1 0.1232235 0.784221
## attr(,"Probability")
## [1] 0.95# LY vs. HS
HS_LY <- re[ , "treatHS:treatLY.sire"] /
sqrt(re[ , "treatLY:treatLY.sire"] *
re[ , "treatHS:treatHS.sire"])
plot(HS_LY)
median(HS_LY)## [1] 0.5620322HPDinterval(HS_LY)## lower upper
## var1 0.2035833 0.8322286
## attr(,"Probability")
## [1] 0.95